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Abstract 

We study the effect of a magnetic field on the behaviour of a slender conducting elastic structure, 
motivated by stability problems of electrodynamic space tethers. Both statical (buckling) and 
dynamical (whirling) instability are considered and we also compute post-buckling configurations. 
The equations used are the geometrically exact Kirchhoff equations. Magnetic buckling of a welded 
rod is found to be described by a surprisingly degenerate bifurcation, which is unfolded when both 
transverse anisotropy of the rod and angular velocity are considered. By solving the linearised 
equations about the (quasi-) stationary solutions we find various secondary instabilities. Our 
results are relevant for current designs of electrodynamic space tethers and potentially for future 
applications in nano- and molecular wires. 
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FIG. 1: Experimental setup for a conducting wire. 

I. INTRODUCTION 

A straight current-carrying wire held in tension between pole faces of a magnet is well 
known to buckle into a (roughly) helical configuration at a critical current (see Figure [I]). A 



photograph of this phenomenon is shown in Section 10.4.3 of 2J], where a linear stability 



analysis is carried out for a simple string model. (A string is here meant to be a perfectly 



flexible elastic wire.) The problem was studied by Wolfe [2l| by means of a rigorous bifurca- 
tion analysis for a (nonlinearly elastic) string suspended between fixed supports and placed 
in a uniform magnetic field directed parallel to the undeformed wire. He found that an 
infinite number of solution branches bifurcate from the trivial straight solution, much like in 
the Euler elastica under compressive load. Minimisation of the potential energy indicated 
that the first branch of bifurcating solutions is stable while all other branches are unstable. 



In a subsequent paper Wolfe 22| extends the analysis to a uniformly rotating (whirling) 
string and shows again the existence of bifurcating branches of whirling non-trivial solutions. 
This result was further extended by Healey 7J using equivariant bifurcation theory in order 
to deal with the symmetries of the problem, which caused the bifurcations to be degenerate. 



Wolfe also considered a conducting rod in a uniform magnetic field [23|. In addition to 
extension a rod can undergo flexure, torsion and shear, and for the case of welded boundary 
conditions it was found that in certain cases bifurcation occurs, with the usual infinity of 
non-trivial equilibrium states. All the works cited above were content with showing the 
existence of bifurcating solutions and did not study their post-buckling behaviour. The 
Hamiltonian structure of the equations for a rod in a magnetic field was investigated in 17 1 
where it was shown that in the case of an isotropic, inextensible and unshearable rod the 
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equations are completely integrable. 

The study of strings and rods in a magnetic field is of great interest to space tethers. 
Although space tethers in the last 20 years or so have become a well established concept 
in astrodynamics 2I, new designs continue to be proposed that hold great potential for 
future space applications. A space tether is a long cable used to connect spacecraft to 
each other or to other orbiting bodies such as space stations, boosters, payloads, etc. in 
order to transfer energy and momentum thus providing space propulsion without consuming 
propellant. These tethers have been studied as elastic strings (e.g., [lO|) and as dumb- 
bell systems (e.g., [25|]). An important class of space tethers is formed by the so-called 
electro dynamic tethers (ETs). These employ the earth's magnetic field and ionospheric 
plasma to generate a current, according to Faraday's Law, and hence thrust or drag forces 



without expending chemical fuel. An exainp 



prototype of the European Space Agency 
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e is the Short Electrodynamic Tether (SET) 



191 ]. This tether system, which comprises a 



central module from which two tethers each about a hundred metres long extend, is designed 
to operate at high inclination and in low orbit. Due to the shape of the earth's magnetic field 
this means that the desirable orientation of the tether is the horizontal one (i.e., with the axis 
of smallest moment of inertia normal to the orbit plane). The gravity gradient ordinarily 
causes a tether to drift to the stable radial position. Therefore, in order to keep the system 
in the horizontal position, an axial spin velocity is applied for gyroscopic stability [sl. [isl. [l8|. 
This requires significant torsional and bending rigidity of the tether which therefore has to 
be modelled as a rod, not a string. The applied spin causes large deformations that present 
stability issues similar to those in an unbalanced rotor system, which have been studied in 
previous work |l^, Q . 

In this paper we apply large- deformation rod theory to study the effect of a magnetic field 
on the behaviour of a slender conducting elastic structure possibly subject to end forces. 
Both statical (buckling) and dynamical (whirling) instability are considered and we also 
cornpute post-buckling configurations. The work extends the stability analysis of the SET 
in 19| by including the effect of the magnetic field. This effect was considered small for the 
SET but that need not be true for other tether designs. For instance, very long and flexible 
tethers subject to boundary conditions that are not too restraining (e.g., no big end masses) 
might well operate in the region of the first magnetic buckling instability. 

We consider welded boundary conditions. These are appropriate for tethers with suffi- 
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ciently large attached end masses (see Figure [2]). Unlike in string buckling the rod does 
not require a tensile force in the trivial state, but we allow for such an applied force as 
well. The pertinent dimensionless parameter that governs buckling measures the product of 
current and magnetic field against the bending force. We find that magnetic buckling of the 
welded rod is described by a remarkably degenerate pitchfork bifurcation. Wolfe considered 



welded boundary conditions in his statical study in [23|] and encountered degeneracies (even- 
dimensional eigenspaces) because of rotational symmetry of the problem, but we show that 
further complications occur, involving branches connecting the bifurcating branches. 

We also study steady whirling solutions for which we introduce a rotating coordinate 
system. This extends Wolfe's analysis of whirling strings to whirling rods. It is found that 
applied spin resolves the degeneracies of the pitchfork bifurcations, provided that the rod 
is transversely anisotropic. We perform a stability analysis by computing eigenvalues of 
the linearised boudary-value problem about a (quasi-) stationary solution. For this we use a 
continuation (or homotopy) approach that takes advantage of the fact that exact expressions 
for the (imaginary) eigenvalues can be obtained in an appropriate limit (no spin, no magnetic 
field). The eigenvalues in this limit are then traced as system parameters are varied. 

The paper is organised as follows. First, in Section 2, we give more details about the 
tether application and discuss the effect of the earth's magnetic field on an electrodynamic 
tether. Then we present our rod mechanics formulation in Section 3, in which the magnetic 
field enters the force balance equation through the Lorentz body force. For the study of 
whirling solutions the equilibrium equations are transformed to a coordinate system rotating 
at constant angular velocity (Section 4). The linearisation is presented in Section 5 and our 
continuation approach to stability analysis is discussed in Section 6. Results are presented in 
the form of bifurcation diagrams in Section 7, both for the statical and dynamical case. We 
also pick up some secondary bifurcations with associated loss of stability. Finally, in Section 
8 we comment on the implications of our results for electrodynamic tethers and draw some 
conclusions. 

II. ELECTRODYNAMIC TETHERS AND THE EARTH'S MAGNETIC FIELD 

Electrodynamic tethers are electrical conductors that interact with the geomagnetic field 
in such a way that an electromotive force (e.m.f.) is generated along the tether due to 



4 




FIG. 2: An electrodynamic tether in the Earth's magnetic field. 



Faraday's Law [1V\. The electrical circuit is closed by means of two contactors attached to 
the ends of the tether which interact with the surrounding ionospheric plasma and allow a 
current to flow. 

Figure [2] shows an electrodynamic tether connecting two satellites. Let {i,j,k} be the 
orbital frame, which may be assumed to be inertial. This is due to the fact that external 
forces on the tether system, that is, electrodynamic forces, are small enough to consider 
a circular Keplerian orbit which allows to decouple translational orbital dynamics from 
structural deformation dynamics [ll]. Besides, since the angular velocity of the tether 
around the earth is much smaller than the typical spin velocity about its axis, rotational 
orbital dynamics can also be decoupled from structural dynamics of the tether subjected to 
spin rotation, which is the main issue of the present paper 19|. In such a case, the tether 
is assumed to travel with a velocity v in the direction of i. The e.m.f. between the ends 
induced by this motion is given by 

E = Jivx Bo) ■ dl, (1) 

where Bq is the magnetic field and dl is a differential along the length of the tether. Because 
the tether is part of a closed circuit a current / will flow in the direction of increasing E 
and the system functions as a generator. This current in turn gives rise to a Lorentz force 
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Fl through 

dFL = dlx Bo. (2) 

This force can be used to drag the system without expending chemical fuel {2!, . Alter- 
natively, if a current is forced against the e.m.f. the system becomes a motor boosting itself 
to a higher orbit. 

If we denote the position co-ordinates of the tether's central axis relative to {i,j, k} by 
(X, Y, Z), then the current vector, which is directed along the tangent of the tether, can be 
expressed as 

The maximum force is generated when / and Bq are perpendicular. In the ET opera- 
tion conditions both vectors will in general not be perpendicular because the tether is not 
perfectly straight and the magnetic field lines will not be perpendicular to the tether over 
its entire length. In order to represent this imperfection let us therefore assume that the 
magnetic field has an extra component in the k direction, 

Bo = {0,B^,B2f, (4) 

where the desired component of the field Bi ^ i?2. Introducing Q and (jl]) into ([2]), the 
differential of the Lorentz force is found to be 

The term —^Bilds in the i component opposes the motion and drags the ET system, as 
intended. However, the crossed B2 terms in the i and j components (due to the imperfection 
B2) wi 
flights 2|,bJ. 

In this paper we study the interaction of elastic and electromagnetic forces in a conducting 
rod and the tendency to generate three-dimensional coiled configurations (for instance as a 
result of buckling at critical loads). 



Intend to coil the tether. This undesirable effect has been reported in some tether 



III. THE ROD MECHANICS MODEL 



We describe the elastic behaviour of a conducting cable by the Kirchhoff equations for 
the dynamics of thin rods. The rod is assumed to be uniform, inextensible, unshearable 
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FIG. 3: Cosserat model of a rod. 



and intrinsically straight and prismatic. The inextensible and unshearable assumptions are 
appropriate for thin rods with low external forces (electrodynamic), as the tether. For the 



background of the Kirchhoff equations the reader is referred to [1|, 13|. These equations were 
also used in [19] to analyse the dynamics of the SET. To these equations will be added the 
Lorentz force to account for the electromagnetic interaction. 

Let X denote the position of the rod's centreline and let {di, d2., d^} be a right-handed 
orthonormal frame of directors (the Cosserat triad) defined at each point along the centreline. 
Since the centreline is assumed to be inextensible we can take ds in the direction of the local 
tangent: 

x\s,t)=d^{s,t), (6) 

where the prime denotes differentiation with respect to arclength s measured along the 
centreline, and t is time. The directors di and d2 will be taken to point along the principal 
bending axes of the cross-section (see Figure [3]). The unstressed rod is taken to lie along 
the basis vector A; of a fixed inertial frame {*, j, k}. 

Looking at Figure [3] we note that the position vector of an arbitrary point of the rod can 
be expressed as 

X(s, 6, 6, t) =x{s, t) + eidi(s, t) + 6^2(5, t) 

(7) 

=a;(s,t) + r(s,^i,6,t), 

where (.^1, ^2) are the components of r in the cross-section relative to {(ii(s), d2{s)}. The rod 
is thus viewed as a set of infinitesimal slices centred at all s. A one-dimensional description 
will be obtained by averaging of forces and moments over each cross-section. The internal 
traction, which is the projection of the stress tensor onto the cross-sectional plane, is given 
by a force which we denote by / = f{s,^i,^2,t) (see Figure [3]). The resultant elastic force 



exerted in a section S{s) is given by 

F{s,t)= [ f{s,^,,^2,t)dS, (8) 

Js{s) 

where dS is an infinitesimal area element. This force can be expressed in the director basis 
as -F = Yl^=i ^i- The resultant moment in the section S{s) is given by 

M{s,t) = [ r(s,6,6,t)x/(s,6,6,t)d^, (9) 

Js(s) 

and will be expressed as M = Yl^=i Midi- 

The rod is assumed to carry an electric current for which we can write 

I = Ix' = Ids. (10) 

Here we have assumed the current to have the same direction as the rod, which is consistent 
with a one- dimensional rod theory. It amounts to the assumption that the cross-section 
of the conducting wire is small enough to make currents within the cross-section (eddy 
currents) induced by the motion negligible. The current / interacts with the magnetic field 
-Bo to generate a (Lorentz) body force given by 

Fl = Ids X -Bo. (11) 



2l| we assume the magnetic field to be uniform and directed along the unstressed 



Following 
rod, i.e., 

Bo = Bok, (12) 

modelling the undesired component of the magnetic field, B2 in equation (jl]). 

The balancing of forces and moments across an infinitesimal rod element then yields 
the following set of partial differential equations: 

F' + IBods xk = pAx, (13) 

M' + dsx F = pihdi X di + /ids x ds), (14) 

where p is the (volumetric) mass density, A the cross-sectional area, Ji and I2 the second mo- 
ment of area of the cross-section about di and ds respectively, and ( ) denotes differentiation 
with respect to time. 
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For a closed system of equations these balance equations need to be supplemented by 
constitutive relations that characterise the material behaviour of the rod. We assume the 
rod to be made of homogeneous isotropic linear viscoelastic material so that stress-strain 
relations, based on a model by Valverde et al. [19;], are 

Ml = E/i(/€i + 7,«:i), 

M2 = Eh{K2+l.k2), (15) 
M3 = G'J(/t3+7»/^3), 

where ki and K2 are the curvatures about di and <i2, respectively, while K3 is the twist about 
ds. The constant 7i, is the viscoelastic coefficient of the material, E is Young's modulus, G 
is the shear modulus and J is the second moment of area of the cross-section about 0^3. We 
shall assume that the section is symmetric with respect to the principal axes, in which case 
J = h + h. 

The are the components of the curvature vector 

3 

K = ^ di, (16) 

i=l 

which governs the evolution in space of the frame of directors as one moves along the 
centreline: 

d[ = K,xdi (i = 1,2,3). (17) 

The constitutive relations f[T^ can be used to replace the Ki in f[T7l) by moments, after 
which the equations ([6]), ( fT3l) . IHM and (|T71) form a system of 18 differential equations for 
the 18 unknowns {x, F, M, di,d2, d^). 

Remark: In general when a conducting wire moves in a magnetic field an additional elec- 
tromagnetic induction effect occurs which opposes the motion. The electromotive force as 

CI 

a result of this effect is proportional to the rate of change of the enclosed magnetic flux [9| . 
However, in the case of a steadily whirling wire in a uniform magnetic fleld the enclosed 
magnetic flux does not change and the effect is zero. 

A. Equations of motion in the moving frame 

Since we are interested in steadily rotating solutions we transform the equilibrium equa- 
tions (|T3|) and (|T4l) to a co-ordinate frame {61,62,63} that rotates with constant angular 
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velocity u) = uok about the k axis (and the axis of the rod in its trivial unstressed state). 
Noting that the derivative with respect to time of an arbitrary vector V{s,t) is given by 



dV{s,t) 



dt 



dV{s,t) 



dt 



+ a;xV(s,t), (18) 



where 4: 1 indicates the derivative with respect to time in the inertial frame and 4i I stands 

dt\i ^ dt\m 

for the derivative with respect to time in the moving frame, the equations f[T^ and f[T^ 
expressed relative to {61,62,63} become 

F' + IB^d^ X 63 = pA{x + 2a;xa3 + cjx(ci;x x)), (19) 
M' + ^3 X F = phidi X di + 2di x (uj x di) + iuj ■ dAidi x uj)) 

(20) 

+ p/i(d2 X ^2 + 2^2 X (cj X ^2) + (<^ ■ d2){d2 X a;)). 
Steadily rotating (whirling) solutions satisfy the equations f|T9|) and fl20l) with the dotted 
variables set to zero: 

F' + IB^d^ X 63 = pAcj X (cj X a;), (21) 

M' + ds X F = p/2(u; • <ii)(di x cj) + p/i(cj ■ d2)(d2 x uj). (22) 

The other equations ([6]) and (ITTI) do not change their form, but all vectors are now to be 
considered as expressed relative to the rotating frame {61,62,63}. Statical solutions are 
simply obtained by setting uj equal to zero. 

For a well-posed problem the final 18 ODEs require 18 boundary conditions to be speci- 
fied. 



B. The boundary conditions 



We follow [23| and consider welded boundary conditions. These conditions also describe 
an electrodynamic tether that is welded to the end contactors or modules if these bodies are 
sufficiently massive (see Figure [2]). We assume the rod to be fixed at s = L and to be able to 
slide along 63 = A; at s = where a controlled force T is applied (positive for tension). L is 
the length of the rod. Writing x = xci + ye2 + ze^ a consistent set of boundary conditions 
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is thus given by 



at s = and 



x{0,t) = 0, 

?/(o,t) = o, 

F(0,t)-e3 = T, 
d3{0,t)-ei^0, 
d3(0,i)-e2 = 0, 
di(0,t)-e2 = 0, 



x(L,t) = 0, 
l/(L,t) = 0, 
z{L,t)^L, 
d3(L,t)-ei^0, 

dsmt)- €2 = 0, 
di(L,i) -62-0, 



(23) 
(24) 
(25) 
(26) 
(27) 
(28) 



(29) 
(30) 
(31) 
(32) 
(33) 
(34) 



at s = L. To these conditions we have to add conditions that ensure the orthonormahty of 
the director basis, for which we can take 



di{0,t)' 


■di(0,t) 


= 1, 


d2{0,t)' 


■d2{0,t) 


= 1, 


ds{0,t) 


■ds{0,t) 


= 1, 


di(0,t)' 


■ d2{0,t) 


= 0, 


di{0,t)' 


■ds{0,t) 


= 0, 


d2{0,t) 


■d3i0,t) 


= 0, 



(35) 



for a total of 18 boundary conditions, as required. Note that these conditions imply that at 
s = and s = L the director frame {di, d2, d^} is aligned with {ei, 62, 63}. 
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IV. NONDIMENSIONALISATION 



We make the system of equations dimensionless by scaling the variables in the following 

way 



F = F 




t = too. 



T = T- 



_ S _ IV X 

c, s = -e[0,i\, uj = —, x = -, 

Li (jjf. L/ 



(36) 



M = M- 



n = hL. 



fEh' ' ' fEh' - -fEh 
Here Uc is a reference bending natural frequency which, through tuning of the numerical 
constant /, can be adapted to the particular boundary conditions at hand and the natural 
mode considered. In the analysis in Section |VT] we shall take / = 1, but since in the welded 
case we consider realistic data for electrodynamic tethers we shall take / = 500.5639 in that 
case so as to get convenient numbers when presenting our numerical results. This value for 
/ corresponds to the first bending natural frequency of a welded- welded beam about di. 

With this nondimensionalisation the equations become (dropping the overbars for sim- 
plicity and letting a prime denote ^): 



F' + Bd^ xk = x + 2(jjxx + (jjx{(jjxx), 



M' + d',x F = P 



R{di X di + 2di X (jj X di + {(jj ■ di){di x cj)) 



+ (^2 X d'a + 2d2 X Lj X d2 + {uj ■ d2){d2 x uj)) 
x' = ds, 



d', = K, X di 



(37) 

(38) 
(39) 
(40) 



and the constitutive relations can be written as 



M 



1 
7 



r(i + R) 

[ki + 7Ki)(ii + R{k2 + 7/^2)^2 H ^ (/^s + 7/^3)^3 



where the dimensionless parameters are 

h 



^ Ji' fEh' E' 



7 = Iv^^c, 



(41) 



(42) 



and (y — 1) is equal to Poisson's ratio. For the boundary conditions we can still use ( !23|) 
to fl35l) if we assume that they now refer to dimensionless variables and that the right-hand 
conditions are imposed at s = 1. 
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V. PERTURBATION SCHEME 



We consider whirling solutions (relative equilibria) that are stationary in the moving 
frame {ci, 62, 63}. Such solutions are found by solving the set of equations flTrj) - fHT]) with 
the dotted variables set to zero (thus obtaining an ODE). To study their stability we linearise 



the full PDE (!37|) - pT]) about these whirling solutions. We follow the approach in [19|, which 
is similar to approaches in s], ^. The stability of statical (non- whirling) solutions can be 
investigated by simply setting the angular velocity u to zero. 
We start our perturbation analysis by writing 



d,{s,t) = d%s) + 5d\{s,t) + 0{5^ 



1,2,3, 



(43) 



where <i^(s) are the basis vectors of a quasi-stationary solution, d*(s, t) are the basis vectors 
of a time-dependent perturbation and 5 is a small bookkeeping parameter introduced to 
separate scales. Note that, in order to preserve orthonormality to 0{5) {di-dj = 5jj -|-0(5^)), 
we must have 



3 

E 



Ms,t)d%s 



1,2,3, 



(44) 



where the matrix Aij is skew-symmetric and can be written as 



\ 



v 



^3 —a2 
-^3 ai . (45) 

a2 —ot\ 

Thus, the nine components of the director basis perturbation are described by only three 
independent parameters, and if we introduce 



cx 



(46) 



(with respect to the unperturbed director basis) then the perturbed director basis can be 
expressed as 

di(s,t) = d°(s) + (5a(s,t) X d°(s) + 0((52), 2 = 1,2,3. (47) 

Using (147|) . the perturbation of an arbitrary vector V = X]i*=i can be written on the 
basis {d^, dg, cZg} as 

3 

V = + + 0(5^) = + + (« X + 0(5'), (48) 



2=1 
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where ()j denotes the component along d- and time and space dependence of the variables 
have been suppressed for the sake of simplicity. 

Applying this perturbation scheme to the PDEs ( l37l) -(l4T i) and the boundary conditions, 
we arrive at an 0(1) nonlinear ODE for the quasi-stationary solutions and an 0{6) linear 
PDE governing their stability. 

A. The 0(1) equations — quasi-stationary whirl 

The 0(1) equations are time- independent. Recalling that a; = uje^, we find the 0(1) 
terms of the linear momentum equation (I57|) . projected on the director basis {d^, dj, dg}, 
to give 



{F^y - F^K^^ + F^kI + 5(^32^21 ~ '^3l'^22 



) = -c.^(x°rf?,+«2) 



,OjO 



{F^y - F^4 + F',n\ 



N 

32^5 



(49) 
(50) 
(51) 



where subscripts are used to indicate components relative to the basis vectors {d^jdgjdg} 
(but the d° components are relative to {ei, 62, 63}). Similarly, the 0(1) term of the angular 
momentum equation fl38l) . projected on the director basis {d'j',d2,d3} gives 



(M?)' = - + + Pu'dU4A, - dU%), (52) 

(M°)' = -^^0^ + /M°M° - + PRu;'dl,{dU% - d?,^"^), (53) 
(M°)' = - fM'.Ml + PRu^dUdUli - dlA2) + Pu:'dl,{dlAi - ^2i^32).(54) 



The 0(1) term of equation (|39|1 can be expressed as 

{x'y = 4, (55) 

and the twist equation ( l40l) by 

(d°)' = /c°xd°, ^ = 1,2,3, (56) 



14 



where = Yl^j=i ^j'^i- "^^^ term of constitutive relations (HT!) can be expressed as 



M° = i< M° = |4 M° = Ii^4 (57) 

which can be used to express the in fl56l) in terms of the moments M^. 

Proceeding in the same way, the 0(1) part of the boundary conditions is given by 

x°(0) = 0, (58) 

/(O) = 0, (59) 

F%0) ■ 63 = T, (60) 

4(0) ■ ei = 0, (61) 

4(0) • 62 = 0, (62) 

d?(0) ■ 62 = 0, (63) 



x°(l) = 0, (64) 

y%l) = 0, (65) 

^°(1) = 1, (66) 

4(1) • 61 = 0, (67) 

4(1) • 62 = 0, (68) 

<(1) ■ 62 = 0, (69) 



4{o) ■ 4(0) = 1, 

4(0) ■ 4(0) = 1, 

4(0) ■ 4(0) = 1, 

4(0) ■ 4(0) = 0, 

4(0) ■ 4(0) = 0, 

4(0) ■ 4(0) = 0, 
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B. The 0{5) equations — linearisation 

The 0{S) part of the hnear momentum equation (137|) can be written as 

{F\s, t))'+Bi{s)F\s, t) + B2{s)x\s, t) + B3{s)a'{s, t) + Bi{s)a{s, t) 

(71) 

=B,{s)x\s,t) + Be{s)x\s,t), 

where the 3x3 matrices Bi{s) are given in the Appendix. Here we have expressed 
relative to {d'j', d^, dg} and a;* relative to {61,62,63}. For the 0{6) part of the angular 
momentum equation (1551) we can write 



(M*(s,t))'+Ci(s)M*(s,t) + C2{s)a'{s,t) + C3{s)a{s,t) + C4(s)F*(s,t) 

(72) 

=C5is)a{s,t) + C6{s)a{s,t), 

where the matrices Ci{s) are again given in the Appendix. is expressed relative to 
dg, dg}. The 9 twist equations (HUjl at 0{5) are reduced to only 3 independent equations 
that relate k* and a as 

K*(s, t) = a'{s, t) + K°(s) X ck(s, t). (73) 
Introducing these relations into the 0{6) part of the constitutive relations gives 

M\s, t) + D,{s)cx'{s, t) + D2{s)cx{s, t) = D^{s)di{s, t) + Di{s)-^^{cx{s, t)), (74) 



where the matrices -Dj(s) are given in the Appendix. Finally, the 0{5) part of equation (!39ll 
yields 

{x\s,t)y = a{s,t)x4{s). (75) 
Applying the perturbation scheme to the boundary conditions at 0{6), we obtain 

(76) 
(77) 
(78) 
(79) 
(80) 
(81) 





= 0, 




= 0, 


F*(0,t)-63 


= 0, 


a(0,t) 


= 0, 


x\l,t) 


= 0, 


a{l,t) 


= 0. 
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After elimination of the k° by means of (157|) . the set of 12 equations (ITTi) . (172|) . (17^ and 
f l75|) together with the 12 boundary conditions fl76l) - fl8Tl) . with appropriate initial conditions, 
form a well-posed initial-boundary-value problem. 

VI. STABILITY ANALYSIS 

Since we are interested in stability of solutions we look for solutions of the 0{6) equations 
of the form 

x\s,t) =x\s)e^\ (82) 
a{s,t) = a{s)e^\ (83) 
F\s,t) = F\s)e^\ (84) 
M*(s,t) = M (s)e^*. (85) 

When these expressions are inserted into (17T i) -(l75 l) a linear eigenvalue problem for a 12- 
dimensional ODE is obtained in terms of the variables {x^,cx.^,F ,M ). The eigenvalue 
A measures the growth of small perturbations and is to be found as past of the solution. 
Eigenvalues come as complex conjugate pairs. A whirling state is unstable if at least one of 
the (in general infinitely many) A's has positive real part. 

To solve a real system of equations we split the eigenvalues and variables (eigenfunctions) 
into real and imaginary parts, X = \r + iXi, x = x^. + ix^, a. = + ia^, F = + iF^ 
and M = + iMi- The equations (17T1)-( 175|) . along with the boundary conditions ( 176|l ~ 
flHTj) . are similarly split into real and imaginary parts. Thus we end up with a doubled 
24-dimensional linearised boundary-value problem. 

A. Stability of the straight rod — magnetic buckling 

The trivial solution of the 0(1) equations (H9|) - fl57|) . representing a straight and untwisted 
rod, is given by 

x{s) = ses, F{s) = -Tes, M{s) = 0, d,(s) = (z = 1, 2, 3), s G [0, 1]. (86) 
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For the statical case {u = 0) without end force (T = 0) the 0{S) equations (17T|) . (1721) . (17il) . 
rSj) about this trivial solution, on inserting fl82l) - fl85l) . can be written as 



X 



A2 B , 



yPx" + —x-—y' = Q 



R R- 

y"" - X^Py" + X^y + Bx' = 0, (87) 

2 \2 p 

M^' - -j^Ms = 0, 

with boundary conditions 

x{0) = x(l) = x'(0) = x'(l) = 2/(0) = y(l) = y'(0) = y\l) = M^(0) = M^(l) = 0, (88) 

while F3 = 0, z = 0. Note that the torsional (M3) modes decouple from the bending {x,y) 
modes. 

To find the statical magnetic buckling loads we set A = 0. The bending equations then 
reduce to 

z""" + ^z = 0, for z = x', (89) 
R 

subject to 

x{0) = x{l) = x'(0) = x'(l) = x""{0) = x""{l) = 0. (90) 

On setting z = e^^'^ we obtain the characteristic equation —k^ + B^/R = with solutions 
fci,2 = ±/3, A:3,4,5,6 = ±/9 (1 ± iV^) /2, where /3 = B^'^/R^'^. Application of the boundary 
conditions fl90l) to the general solution z{s) = Yl^j=i o^jC^^^^ leads to the following critical 
condition: 

x(/9) := 2 cos (3 + cos 2/5 — 2 ^cos ^ + cos cosh — ^ + (2 — cos /5) cosh V^(3 

— V^sin/? sinh a/3/5 — 2v^ ^sin ^ — sin sinh = 0. (91) 

The critical loads correspond to pitchfork bifurcations where non-trivial solutions bifurcate 
from the trivial straight solution. A plot of x (Figure H]) shows that in the welded case the 
pitchfork bifurcations are (doubly) degenerate, as was also found by Wolfe {23I. We stress 
that the above calculation is only possible for the statics case. If 7^ then the x and y 
equations do not decouple and no simple characteristic equation is obtained. 
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FIG. 4: Degenerate zeroes of the characteristic equation for the welded rod. The first three critical 
[i values are 7.332130, 13.613561 and 19.896753. 

B. Eigenvalues for the unperturbed problem (T = 0, 7 = 0, = 0, i? = 0) 

We shall call the case where T = 0, 7 = 0, ct; = and 5 = the unperturbed problem. 
For this problem explicit expressions can be obtained for the eigenvalues of the linearisation 
about the straight solution. The x and y equations in ( IHTl) decouple into two fourth-order 
beam equations: 

R ' (92) 



a;"" _ A^Px" + — X = 0, 



y"" - x'Py" + x'y = o, 



subject to boundary conditions (jSHl) . Since we anticipate imaginary eigenvalues we set 
X = ifj,, X = e^'"^, y = e*'^'*, and find for the x equation 

1/2 



/I 1 \ ' 

/I , 1 \l/2 

h,^ = ±1 l-^fi^P^ + 4fiyR--fi'p] =: ±ib. 



while for the y equation 

1/2 



/I 1 
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where a, b, a, (3 are non-negative real numbers. The general solutions are 



x{s) = Ax sin as + cos as + Cx sinh bs + Dx cosh 6s, 
y{s) = Ay sin as + By cos as + Cy sinh /5s + cosh /5s. 

Application of the boundary conditions fl88p leads to 

(a^ — 6^) sin a sinh 6 = 2a6(l — cos a cosh 6), 
(a^ — /?^) sin a sinh/? = 2a/?(l — cos a cosh/?). 

In order to obtain the eigenvalues these two transcendental equations can be solved nu- 
merically using a Newton-Raphson scheme. Meanwhile, the torsional eigenvalues of the M3 
equation in (157|) are given by 




These are all the eigenvalues for the unperturbed problem. They will be used as starting 
values in the numerical procedure described next. 



C. Numerical procedure 

The main idea is to use the known eigenvalues in the unperturbed problem as starting 
values in a continuation procedure in order to compute the eigenvalues and corresponding 
eigenfunctions for general values of the parameters T, 7, u and B. For this we use the well- 
tested code AUTO (specifically AUTO2000). AUTO solves boundary- value problems 
by means of orthogonal collocation. It requires a starting solution and can then trace out 
solution curves as a parameter of the problem is varied. Bifurcations are detected where 
branches of solutions intersect. At such points AUTO is able to switch branches and compute 
curves of bifurcating solutions. 

Our procedure takes advantage of the fact that A appears only quadratically in the 
linearisation ( jTTl) . ( ff2l) . (ff4l) and (ff5l) if 7 = and u = 0. To explain the method consider 
the typical 0{6) equation 

z"" - X^f{s)z" + X^g{s)z = 0, (93) 
where / and g are functions of the 0(1) solution. Writing z = x + iy, A = + iXi, we can 
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decompose the z equation into 

x"" - (A^ - \\)f{s)x" + 1\\f{s)y" + (A^ - \^^g{s)x - 2X,Kg{s)y = 0, 

(94) 

y"" - {Xl - X^)f{s)y" - 2\XJ{s)x" + (A^ - >^)g{s)y + 2K\g{s)x = 0. 

The important thing to note here is that these equations decouple into two identical equa- 
tions if the eigenvalue is either imaginary (A^ = 0) or real (Aj = 0). 

This suggests the following sequence of steps, involving boundary-value problems of in- 
creasing dimension, to compute eigenvalues of statical or uniformly whirling solutions. 

1. Consider the unperturbed problem of Section IVl Bl and, noting that all eigenvalues are 
purely imaginary, solve the 30-dimensional system of 18 0(1) equations and one 12- 
dimensional system for the imaginary part of the 0{6) equations (cf. the y equation in 
(19^ ). Set Ar- = and use Aj as the continuation parameter in AUTO in order to com- 
pute the eigenvalues (instead of solving the transcendental equations in Section IVIBp . 
These eigenvalues will show up as branching points (BP), or pitchfork bifurcations, as 
eigenvalues by definition are those values for which non-zero BVP solutions exist. By 
symmetry it is only necessary to consider Aj > 0. 

2. Keeping the same 30-dimensional system, switch branches at a BP to compute ('grow') 
the corresponding (imaginary) eigenfunction. Since the equations are linear the value 
of Aj will not change in this run. For later use we monitor the non-zero solu- 
tion by means of some measure ||.||j (not necesarily a proper norm) on the space 
{a;-, CK*, Mi} of imaginary linearised variables. 

3. Now consider the full system of 42 equations (18 0(1) equations and two sets of 12- 
dimensional 0{6) equations (cf. fl94p ) . Fix the measure on the imaginary part 
and release A^ instead in order to compute the real eigenfunction (since the imaginary 
part of the solution is fixed there is only one branch of solutions through the starting 
point and there is nowhere else to go for the continuation but to 'grow' the real 
eigenfunction). Again we monitor this function by means of a suitable measure ||.||t-- 
In this run neither A,, nor A, will change. 

This approach works because the solution obtained in step 2 also solves the full 42- 
dimensional system when the extra 12 variables (i;* , a* , F^, M^) set to zero. This 
is a consequence of the fact that the real and imaginary parts of the 0{6) equations 
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decouple if A,. = 0, as a result of the quadratic dependence of the eigenvalue problem 
on A (cf. (1941)). 

Steps 2 and 3 can be performed for as many of the BPs computed in step 1 as required and 
will give the corrsponding eigenvalues and eigenfunctions. Once these have been obtained 
both measures and ||.||j can be fixed and an extra system parameter such as -B or a; 
released in order to trace the eigenvalues (and hence monitor stability changes) as system 
parameters are varied. (Note that fixing \\.\\r and ||.||j makes sense as eigenfunctions are 
only defined up to a multiplicative factor.) 

The above 3-step procedure is not limited to linearisations about the trivial straight 
solution. It can be applied to any starting solution that has no eigenvalue with both A^ and Aj 
non-zero, as these would not be picked up in step 1. (It is of course no problem if eigenvalues 
become fully complex (e.g., in a Hopf bifurcation) in the course of further continuations.) 
For instance, we find that at the first critical B, given by (!9T|) . the lowest conjugate pair of 
eigenvalues ±Aj goes to zero and becomes a real pair of eigenvalues, signalling a stability 
change of the straight rod. The (first-mode) solution bifurcating at this point is stable 
with all eigenvalues being imaginary and the above procedure can be applied to find the 
eigenvalues. 

We end this section with a few comments: 

(i) There are infinitely many eigenvalues and the above procedure only finds the lowest 
order ones. This is of course a limitation of any numerical scheme. We find that 
eigenvalues vary slowly with system parameters, suggesting that stability is governed 
by the lowest-order eigenvalues. We typically consider 5 or 6 eigenvalues. 

(ii) Note that in steps 1 and 2 above we could not have taken the full 42-dimensional system 
of equations as that would have made the branching points (pitchfork bifurcations) 
degenerate and AUTO would not detect a BP. This is because if A,. = (or Aj = 0) 
the two sets of 12- dimensional linearised equations are identical (cf. flM|) ). 
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TABLE I: Dimensional and dimensionless parameters for the SET. 



L 


100 m 


P 


0.001 


A 


2.879 X 10-11 m^ 


R 


0.5526 


E 


1.32 X 10^1 N/m^ 


V 


0.76923 


Eh 


38 Nm^ 


f 


500.5639 


Eh 


21 Nm2 







VII. NUMERICAL RESULTS 
A. The statics case, w = 

Figures O and [H] show the bifurcation diagram obtained when the magnetic field parameter 
B is varied, using different measures along the vertical axes. The dimensionless parameters 
taken are those of the SET and are listed in Tabled along with the dimensional parameters. 
We observe that at critical values of B pairs of non-trivial solution branches bifurcate, 
thus confirming the degenerate nature of the bifurcation found analytically in Section IVI A[ 
The first three bifurcations occur at B = 0.585 (5P1), 3.747 (5P2) and 11.698 (5P3), in 
agreement with Figure H] when the numerical factor / is taken into account. In addition 
to these primary branches, there are branches, labelled bl2, b34 and b56, connecting these 
pairs. One end of these connecting branches reaches down to the horizontal axis, while the 
other end comes out of a secondary pitchfork bifurcation along the first branch of the pair of 
primary branches. They in fact come in symmetric pairs as is brought out by the measure 
plotted in Figure El The shapes of the rod along the connecting branches form a smooth 
transition between the shapes on the connected branches; see Figure [7] where projections 
onto the {ei-es} and {62-63} planes of bifurcating solutions along branches bl, b2 and bl2 
are shown at constant measurci = 0.04. 

Solutions bifurcating at larger B values have successively more coils. Figure [8] shows the 
three-dimensional shape of a solution on the fifth bifurcating branch. 

Figure [9] shows the evolution of the imaginary and real parts of the first five pairs of 
eigenvalues along the trivial solution, from B = to B > BP2 (i.e., beyond the second 
pitchfork). As the branches are born in pairs at the degenerate pitchfork, two pairs of pure 
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FIG. 5: Bifurcation diagram for the statical case {to = 0). System parameters are those of Table HI 
measurei = |2;'^(s)|ds. 
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FIG. 6: Bifurcation diagram for the statical case (lj = 0) using a different solution measure. 
measure2 = /o x°(s)j/°(s)ds. The primary branches have measure2 = 0. 
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FIG. 7: Projection of the modes corresponding to branches bl, b2 and transition branch bl2 onto 
the {61-63} and {62-63} planes at constant measurei = 0.04. Values of B are: 0.595 (bl), 0.605 
(bl2) and 0.7 (b2). 
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FIG. 8: 3D view of a solution along branch b5. 

imaginary eigenvalues collide at zero and become pure real pairs at i? = 0.585 [BPl). One 
of the eigenvalues of each real pair is positive, meaning that the trivial solution becomes 
unstable a.t B = BPl, as one would expect. A similar further loss of stability occurs at 
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FIG. 9: Evolution of imaginary (a) and real (b) part of the first five pairs of eigenvalues along the 
trivial solution as B is increased from to > BP2. 

B = BP2. 

Figure [10] shows the evolution of the first five pairs of purely imaginary eigenvalues when 
switching from the stable trivial branch to branch bl, which is initially found to be stable. 
The figure reveals that branch bl loses stability at the secondary pitchfork bifurcation at 
B = 1.942, where bl2 connects. At this point the first pair of eigenvalues becomes real, as 
illustrated in Figure [HI This figure also shows that connecting branch bl2 is unstable as it 
too has a pair of real eigenvalues. 

Considering the problem of the ET from a design point of view, it is interesting to 
quantify the value of the geomagnetic field at which the straight tether would buckle. In 
the case at hand, this value is B = 0.587 (see Figure [5]), which in dimensional parameters 
yields IBq = 2.104 x 10"^ — . Noting that the maximum value of the geomagnetic field is 



Bg = 7xlO-^T (lit 
would he I = 1 A 



, and assuming that the maximum current that will flow along the tether 
llj, the maximum expected value for the constant IBg = 7 x 10^^ ^. 



This means that the ET would be below the critical value at any time. Note, though, that 
B goes as the cube of L and for longer tethers, which are common in radially stabilised ETs, 
the critical value may be exceeded resulting in buckling into a coiled shape (cf. Figur e Ml , 
as has been reported in some tether flights such as the PGM and TSS-IR missions 2|, lll |. 
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FIG. 10: Imaginary part evolution of the first five pairs of eigenvalues along the trivial solution 
from i? = to = BPl, then following branch bl beyond the secondary pitchfork. 

B. Whirling solutions (relative equilibria), uj ^ 

We now introduce angular velocity u to the ET and seek relative equilibria, that is, 
solutions that appear static when viewed from the rotating frame {ci, 62, 63}. Figures [T2l[T3] 
and [TH show the effect of increasing angular velocity on the bifurcation diagram of Figure [5l 
The degenerate pitchfork bifurcations get resolved and the connecting branches move up 
along the principal branches. For uj = 0.25 (Figure [T2|) all branches, bl, b2, bifurcate, 
as in the statics case of Figure [5l and the trivial solution is stable up to the first primary 
pitchfork. However, at uj = 0.75 (Figure [T3l) . bl has merged with its symmetric partner at 
-B < and lifted off the horizontal axis, leaving no stable straight solutions. b2 similarly 
lifts off near u = 1.25. Figure [TH shows the bifurcation diagram at u = 2. An extra 
branch of connecting solutions (b22) has appeared, but no further branches have become 
non-bifurcating. 

Figure [15] shows the bifurcation diagram at uj = 2 for the case of an isotropic rod [R = 1). 
Here the bifurcations become degenerate again but the connecting branches do not move 
down. In fact, they turn into vertical branches. We conclude that neither anisotropy (R) 
nor whirl [uj) alone resolves the pitchfork degeneracy; both are required. 

To investigate the stability of non-trivial solutions we follow the eigenvalues computed 
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FIG. 11: Real part evolution of the first pair of eigenvalues along the trivial solution from i? = to 
B = BPl, then following branch bl beyond the secondary pitchfork as well as connecting branch 
bl2. 




FIG. 12: Bifurcation diagram for a whirling anisotropic rod at lij = 0.25. Branches bl and b2 still 
bifurcate (at distinct values of B) and trivial solutions before the first pitchfork are stable. 
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FIG. 13: Bifurcation diagram for a whirling anisotropic rod at cj = 0.75. Branch bl no longer 
bifurcates and all trivial solutions are unstable. 




FIG. 14: Bifurcation diagram for a whirling anisotropic rod at = 2. Branches bl and b2 no 
longer bifurcate. An extra (unstable) branch b22 has appeared. 
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FIG. 15: Bifurcation diagram for a whirling isotropic rod {uj = 2, R = 1) 
(a) (b) 





FIG. 16: Evolution with cu of real (a) and imaginary (b) parts of the first five pairs of eigenvalues 
from branch bl of Figured {uj = 0) with B = 0.5869. 

earlier as uj is increased from 0. Figure [16] shows the evolution of the real and imaginary 
parts of the first five conjugate pairs of eigenvalues when continuing in u from a static 
solution on the first branch (bl) of Figure [5] at fixed B = 0.5869. Note that the first and 
second eigenvalues are almost coincident and can hardly be distinguished in the figure. Both 
eigenvalues acquire positive real parts; hence the static coiled tether becomes unstable as 
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FIG. 17: Evolution with respect to to of the real part of the first three pairs of eigenvalues for a 
solution along bl at 7 = 0.04375 (solid) and 7 = 0.01 (dashed). {B = 1.2.) 

soon as it is being rotated about k. The other eigenvalues move to the negative real half- 
plane as u is increased. The angular velocity is increased to u = 2, which is beyond the 
critical velocity in the B = 0-case where the tether buckles at the first bending natural 
frequency of the beam, 00 = \/R jiol . 

Next we consider the combined effect of both damping and spin on the stability of the 
rod. Figure [T7] and [18] show the evolution of the first three eigenvalues as uj is increased for 
a solution along the first bifurcating branch bl at constant B = 1.2 and two different values 
of the damping coefficient. The value 7 = 0.04375 was estimated by Valverde et al. |l^ 
to be the damping value for a real electrodynamic tether under working conditions. The 
figures show that the solution remains stable at least up to = 4 (higher-order eigenvalues 
do not seem to alter this picture). However, at the lower value of 7 = 0.01 we observe a Hopf 
bifurcation at a; ~ 1.5 where the conjugate pair of eigenvalues (Ai, Ai) acquires a positive 
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FIG. 18: Evolution with respect to to of the imaginary part of the first three pairs of eigenvalues 
for a solution along bl at 7 = 0.04375 (solid) and 7 = 0.01 (dashed). The curves are independent 
of 7. {B = 1.2.) 

real part (at u ^ 2.2 the second pair of eigenvalues (A2, A2) moves into the positive real half 
plane, returning to the negative half plane at u; 6). Also note that damping only affects 
the real parts of the eigenvalues, not the imaginary parts (Figure [T8l) . 

Figure [19] shows curves of these Hopf bifurcations, both of Ai and A2, in the B-u parameter 
plane for various values of 7, as well as the curve of secondary pitchfork bifurcations. The 
latter curve is independent of 7. This makes sense as our (quasi-) statics problem, and hence 
the bifurcation diagrams in Figures [12], [13] and [TH do not depend on 7, and therefore the 
secondary pitchfork bifurcations, where a real eigenvalue goes through zero, must also be 
independent of 7. Interestingly, the Hopf bifurcation curves of Ai are closed, with instability 
inside the curves, while the Hopf bifurcation curves of A2 are open. The closed curves shrink 
to zero at 7 = 0.023, leaving only instability due to Hopf bifurcations of A2 for larger values 
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FIG. 19: Curves of secondary pitchfork bifurcations (dark grey), Hopf bifurcations of Ai (black) 
and Hopf bifurcations of A2 (light grey) at 7 = 0.01, 0.013, 0.016, 0.02 and 0.022 (from outside to 
inside). The stable region is below the dark grey curve and around the black and light grey knees. 

of 7. In the limit 7-^0 the Hopf curves approach the coordinate axes u = 0, B = 0. The 
region of stable bl solutions in the diagram is the connected L-shaped region between the 
secondary pitchfork and Hopf curves. 



VIII. CONCLUSION 



We have shown that whirling current- carrying rods bifurcate under increasing magnetic 
field (or current). For the welded boundary conditions considered here this bifurcation is 
doubly degenerate (even for a transversely isotropic rod). This is because the equations are 
invariant under rotation about the axis (63 axis) of the su ppo rts. This symmetry property 



complicated Wolfe's analysis, both for the whirling string [22| and the statical rod [23]. In 
both models account had to be taken of the variational nature of the problem to prove 
existence of non-trivial bifurcating states. We find that for the rod, eigenvalues become 
simple when rotating states are considered, provided the rod is anisotropic. Neither spin 
nor anisotropy alone resolves the degeneracy. 



Wolfe's analysis [2ll] for non-whirling strings suggests that the first bifurcating branch is 
stable and that all other branches are unstable. We find the same to be true for a welded 
rod, initially, but a secondary pitchfork bifurcation occurs where the post-buckling solution 
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(along bl) loses stability. Unfortunately, our method allows us only to study stationary or 
quasi- stationary (whirling) solutions, so it is not clear what type of solutions occur when the 
(quasi-) stationary solutions lose stability in these secondary bifurcations. To investigate 
this one would have to do simulations based on direct discretisation of the PDEs fITU]) and 

m- 

Associated with the secondary pitchfork bifurcations are connecting orbits (bl2, b34, 
etc.) which interfere with the primary pitchfork bifurcations in the limit uj ^ 0. Another 
interesting bifurcation scenario concerns the merging of symmetrically related branches (at 
B = 0) and the subsequent 'lifting off' of these branches, thereby becoming non-bifurcating. 
In the case of bl (Figure US]) this leads to the disappearance of stable trivial (straight) 
solutions. 

We have performed the, as far as we are aware, first study of the magnetically induced 
post-buckling behaviour of elastic rods. Our results may be of interest to the design of 
static as well as spinning electrodynamic tethers, which often operate in the post-buckling 
regime. As we pointed out in Section 7.1, the critical values of the magnetic field (or current) 
that we find may be close to operating values for certain tether types. Furthermore, axial 
spin tends to destabilise the post-buckling state of the tether through Hopf bifurcations, 
although typical viscoelastic material damping levels go some way towards avoiding Hopf 
instabilities. On the other hand, spin has a stabilising effect in that it tends to push the 
secondary pitchfork bifurcation to higher values of the magnetic field. We have mapped out 
stable regions in the B-u parameter plane (Figure fT9l) . 

We finally like to speculate on another possible application of this work. There is great 
current interest in conducting nanowires. These can either be silicon-based wires, carbon 
nanotubes or metal-coated biological fibres such as proteins, DNA molecules and micro- 
tubules [l6]. In addition there is the ongoing discussion whether or not DNA molecules are 
electrical conductors 121]. All these nanometer-scale structures are believed to have great 
potential as building blocks for future electronic devices. The interaction of these wires with 
magnetic fields could conceivably be exploited to obtain certain desirable properties. 
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Appendix A: Matrices for the linearisation 



The matrices Bi appearing in equation (17T|) are given by 
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Matrices Ci appearing in equation (!72|) are given by 



Ci = B 



1) 



Co 



-M°^ 

-M3O M° 

M° -M° y 



37 



/ r'li r^vi r^i3 \ 

■-^3 '-^3 "-^3 

/^21 n22 nlZ 

^■i 

r^?,\ /^32 ^33 

•-^3 '-'3 '-'3 J 



where 





= Mi4 


+ - P{uj ■ dl){dl xuj-d\)- P{u ■ 


d°)Kxa;.<), 


cf 


= {Mt)' 


- - PR{uj ■ cZ?)(d° xcj-dl), 




cf 


= -(M° 


)' - Ml 4 - Fl + PR{uj ■ dl){dl xuj-dl] 


+ P(Lc;-d?)(d° X 


cf 


= -{M's 


y-MlK\ + P{u.dl){dlxu.dl), 




Cf 


= m^kI 


+ M°K° + PR{u} ■ dl){d\ X a; ■ d°) + PR{ 


■ d\){dl X UJ ■ 


cf 




- M^kI - F° - PR{uj ■ dl){d\ xuj-dl)- 


-P(u;-d°)(d°xa;-cZ?), 


cf 


= 


- MIk\ + F° - P(w • dl){dl X a; ■ 




Cf 


= -{Ml 


)' - MIkI + F° + PR{uj ■ xuj-dl) 




cf 


= 


+ M°K° + P(l - P)(lj ■ d°)(d° xuj-dl)- 


FP(1 -P)(u; ■cZ°)(d° X a; 



0-10 
1 




1 
OP 
(l + R) J 



X X d?) ■ 
P(cZ° X X d°) ■ d° 



^ (d° X X d°) ■ d° 
3 = 2P -P(d° X X d°) ■ d^ 

^ (d^ X X dl) ■ d° -P(d° X X d°) ■ d° P(d? X ^ X d°) ■ d° - (d° X u; X d?) • dj 

Matrices Di appearing in equation fl71|) are given by 
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All the in the above can be expressed in terms of the moments Mf by means 
constitutive relations (E 
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